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I. INTRODUCTION 


Phonon energy dispersion relations (EDR) and the 
density of states (DOS) are a fundamental physical prop¬ 
erty of a solid especially for determining the mechanical, 
thermal and other condensed matter phenomena. Al¬ 
though most of the interest in the past has been in three 
dimensional (3D) systems, for example diamond, with 
the discovery of graphene [T] interest in systems of lower 
dimensionality has arisen, including graphene (2D) [2j|3], 
graphene tubules (ID) [4] and fullerenes (OD) [5]. 

At the same time, the development of superconduct¬ 
ing microwave billiards mm and of photonic crystals 
ElE] combined into ’’microwave photonic crystals”, or 
generally, of ”artificial graphene” pTH2Q] . has opened the 
way for simulations of crystal vibrations in systems of 
lower dimensionality (2D, ID, OD) with a finite number 
of units. 

Calculations of the EDR and the DOS in infinite sys¬ 
tems have been done using a variety of methods, includ¬ 
ing tight-binding force-constant models [H [211 HI] and 
molecular dynamics methods. Here we introduce a new 
method, based on the algebraic theory of molecules [23] . 
which reproduces all previously known results for the fun¬ 
damental vibration v = 1 of crystals, and, in addition, 
can be extended to overtone vibrations. The method is 
particularly well suited for analyzing the EDR and DOS 
of microwave photonic crystals. 

With this method, we calculate the EDR and DOS of 
the fundamental (i; = 1) and overtone {v = 2) transverse 
vibrations of ID and 2D crystals. For 'z; = 1 we also 
rederive the EDR and DOS analytically, determining the 
location of (van Hove) singularities and (Dirac) zeros in 
transverse vibrations of ID crystals and 2D crystals with 
square and hexagonal lattices. The latter is then com¬ 
pared with experimental results in microwave billiards. 
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We are then able to determine the scaling behavior of 
the singularities. 


II. ALGEBRAIC METHOD 
A. One-atom 


We consider a one-dimensional (Id) vibration of an 
atom (for convenience in the ^-direction), and assume 
that it is bound by a potential 


RoA(A-l) 

1 2 ’ 
cosh az 


( 1 ) 


as shown in Fig. The energy levels in this potential 



FIG. 1. The Poschl-Teller potential, Eq. Q, with a = 1, Vb = 
1, A = 3. 

can be described by the simple formula [24] 

E{v) = EoEA{v- , (2) 

where ^ (also written, for integer A, as ^) de¬ 

scribes the anharmonicity, A = 2Vo{X — 1) (also written 
as huj) the oscillator frequency, and Eq = — Vb(A — 1)^ 
the zero-point energy. It was shown in [25] that the 
Poschl-Teller potential Eq. 0 can be associated with 
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the algebra of = u{2). Its eigenstates are described by 
representations of u{2) D 5o(2), labeled by \N^v) where 
= N=eYen or odd) is the vibra¬ 

tional quantum number, and N is the so-called vibron 
number which determines how many bound states are in 
the potential. The Hamiltonian can then be written in 
algebraic form as 


H = E^^ AC, 


( 3 ) 


where Eq is the zero point energy, A a scale and the 
operator C, called Casimir operator, has eigenvalues 


{N,v\C\N,v) 


4 V- — 


( 4 ) 


From Eq. Q, one can see that N controls the anhar- 
monicity of the potential C = For N ^ oc, the spec¬ 
trum is harmonic and has an infinite number of bound 
states, V = 0,l,...,oo. The potential Eq. 0 is therefore 
well suited to describe both harmonic and anharmonic 
vibrations. 


B. Many-atoms 

In the case of many atoms at location i = l,...,n, 
bound by Poschl-Teller potentials, the algebra is ^ = 
ui{2) 0 U2{2) 0 ... 0 Ui{2) 0 ... 0 Un{2) [26l |27]. The 
algebraic Hamiltonian for n uncoupled oscillators is 


H = En 




( 5 ) 


i=l 


where Eq is the zero-point energy and Ai a scale. Its 
eigenvalues are 


E = Eo-4j2AiC-^ 


i=l 


If all the oscillators are equivalent, A^ = A, Ni 
and 


E = Eo-4AY^ 


i=l 


N 


(6) 

N, 

( 7 ) 


In the harmonic limit, Ni ^ oo, the matrix elements 
of the operators Ci are the familiar matrix elements of 
the number operator hi = 6^6^, where bl{bi) are boson 
creation and annihilation operators of a phonon at site i, 
multiplied by —4. 

The phonons (vibrons) interact with a diagonal and 
an off-diagonal interaction. The diagonal interaction is 
written in terms of the normalized Casimir operator Cij 
of the combined S0i(2) 0 S0j(2) algebra as 


(Ei , Vi 5 JVj , Vj I Cij I Ei , Vi 5 JVj , Vj ) 


( 8 ) 
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FIG. 2. Spectrum of equ ival ent harmonic oscillators with di¬ 
agonal interaction, Eq. (11), with Eq = 0, A = —1, N = 
10, A' = -0.1. 


In practical calculations, it is convenient to subtract from 
Cij a contribution that can be absorbed in the Casimir 
operators of the individual modes i and j given in Eq. 
thus considering an operator CF whose matrix elements 
are 


(Ni, Vi; Nj,Vj jC'j I Ni,Vi; Nj,Vj 


= -4 


Ni 


\ - ’Z: ^25 - ’ J 5 / 

m 2 / , n 2 ' 


( 9 ) 


Vj (Vi+VjY 


Nj (Ni + Nj) 


The Hamiltonian for n oscillators diagonally coupled is 

n n 

H = Eo + Y,AiCi + Y, 00. (10) 

i=l i^j 

If all oscillators are equivalent, Ni = N, Ai = A, A'- = 
A', and the eigenvalues are 


E = Eo - 4H 

n 

-«T 


E 

i=l 


N 


v^ 
~ N 


v] {viEVjY 


N 

N 


2N 


( 11 ) 


The diagonal interaction represents the modification in 
the Poschl-Teller potentials at each site due to the other 
sites. It is usually small and it can be neglected. The 
spectrum of states of equivalent oscillators with diagonal 
interaction is shown in Eig. for up to two phonons, 
Vi = 2, both for harmonic {N oc) and anharmonic 
oscillators. The degree of degeneracy of the levels is given 
below in Eq. (19) and is shown to the left of the levels. 


The off-diagonal interaction is written in terms of op¬ 
erators, called Majorana operators M^j, with matrix ele- 
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ments [isiiiniiiT] 

{Ni,Vi + l\Nj,Vj - 1 \Mij\Ni,Vi;Nj,Vj) (12) 

= - [vjivi + im - Vim - vj + 1 )]'/' 

{Ni, Vi - 1; Nj,Vj + 1 \Mij\Ni, vf, Nj,Vj) 

= - [vm + 1) W - v,){Ni -Vi + 1)]^/" /[NiNjm- 
The Major ana operators have also diagonal matrix ele¬ 
ments given by 

{Ni,Vi;Nj,Vj \Mij\Ni,Vi-,Nj,Vj) ( 13 ) 

= {viNj + VjNi - 2viVj) / . 

If the oscillators are equivalent, Ni = their contribu¬ 
tion, (^Vi + Vj — ^ can be reabsorbed into Eq. ( jll| ) 

so it will not be considered further. 

The operator 

n 

M = J2Mij ( 14 ) 


is Hermitian and its matrix elements are real and sym¬ 
metric. It is called Majorana (or exchange) operator 
because it was introduced by Majorana in 1993 in the 
context of nuclear physics [28] . Within the algebraic the¬ 
ory of molecules [23] it was introduced in 1991 [26]. Its 
mathematical definition is that of the invariant operator 
of 1^1 (2) 0 1 ^ 2 ( 2 ) 0 • • • 0 i4n(2). In recent years, operators 
similar to M have been used in the context of condensed 
matter physics. 

The elements Mij of the Majorana operator annihi¬ 
late one quantum of vibration at site i and create one at 
site j, or vice versa. In the harmonic limit, Ni m 00 ^ 
the matrix elements of Mij are the familiar matrix ele¬ 
ments [vj{vi 0 1)]^^^ of the operator blbj, where b] and 
bi are boson creation and annihilation operators of a 
phonon at site i, respectively, multiplied by a minus sign. 
Because of the sum in Eq. (14), M can be written as 

Yl^>i ’ which shows explicitly the Hermitic- 

ity of M. The operator M splits the degeneracy of the 
states in Eig. determines the distribution of eigenval¬ 
ues and thus the level densities for the fundamental vi¬ 
bration V = Vi = 1 and the overtones 'z; = 2, 3,..., and 
generates the normal modes of vibrations. 

The total Hamiltonian for Id vibrations of a system of 
n atoms at locations Vi^i = 1, ...,n, can be written as 


n n n 

F = So + E + E 00 + E 

i=l i^j i^j 

and for n identical atoms as 

n n n 

H = Eo + AY,Ci + A'Y,C'ij + Yl ■ ( 16 ) 

i=l i^j i^j 

This Hamiltonian applies to any system. In sections Ing 
[rv|and[v| we will consider three cases: (1) linear lattice 
(ID), (2) square lattice (2D) and (3) honeycomb lattice 
(2D). 


III. LINEAR LATTICE 


The algebraic theory of linear lattices, Eig. [^ was in¬ 
vestigated in [29]. The most general algebraic Hamilto- 


• - • - ^ - • - • 

1 2 i j n 


FIG. 3. The linear lattice. 


nian for this problem can be written as 
= A ^ C,- + A(') 

3 3 3 

+ (17) 

3 

with j = 1,2, ...,n, and where is the strength of 
the nearest-neighbor interaction, of the next-nearest 

neighbor, of the next-to-next one, etc. 

The Hamiltonian Eq. is of the general form 

H = Aj2Cj + Y.mMjj,. (18) 

3 37 ^ 3 ' 

Denoting by \vj) the number of quanta at each site j, the 
basis for the diagonalization of Eq. ( jlSj ) is: 

-z; = 0 vacuum |0) = |0,0,..., 0) 

v = l n basis states |lj) = |0,0,..., Ij,..., 0) 

v = 2 n basis states |2j) = |0,0,..., 2j,..., 0) 

V = 2 basis states |lj, 1^/) = |0, 0,..., Ij,..., 1^/,..., 0) 


(19) 

The states with v = 1 constitute the fundamental vibra¬ 
tion, the states with v = 2 the overtone and combination 
modes of the linear lattice. The matrix elements of the 
operators Cj are 


{vj \Cj\vj) 



N^oo 


-4vj 


N ’ 


1/2 


N^oo 




( 20 ) 


In the harmonic limit, N m 00 , the algebraic Hamilto¬ 
nian Eq. (18) reduces to the boson Hubbard Hamiltonian 


H = -iAj2 - E (21) 

3 37 ^ 3 ' 


where 6^ creates a quantum at site j. Eor identical oscil¬ 
lators and nearest-neighbor interaction, it can be written 
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as 


H = -iAj2ni-XYl + ^*^ 1 + 1 ) ’ ( 22 ) 


with, in the notation of Ref. [30], /i = 4^4 and t = X. 
Since the operator Cj is diagonal in the basis 


Eq. (19) and the operator ^j^jf Mjjf conserves the total 
numl^ of quanta v = Vj , the diagonalization of H 
splits into blocks. H can then be diagonalized in the sub¬ 
spaces with V = V = 2^ ..., of dimension n, 

It is of interest to note that the algebraic method can 
also describe more complex situations such as those in 
which the algebraic Hamiltonian is 

H = AY^Ci + BY,Cf + B'Y,CiCj + xY,Mij, (23) 






which is equivalent to the boson Hubbard Hamiltonian 
with in-site interactions 


H = hi 


uY^m 


.i{hi-i)+vY 




— t 




(24) 


The diagonalization of this Hamiltonian can be done in 
the same way as before since the addidional terms do not 
modify its block-diagonal form. 


A. Solutions 

1. The fundamental vibration v=l 


If the number of sites is large, one can replace by 

a continuous variable x which runs from 0 to tt and Ek 
by a continuous variable E which runs from a — 2f3 to 
a -1- 2/d. The DOS of the v = 1 vibration is then 


p{E) = 


1 

2/d sin arccos( ) 


(29) 


This result has been verified by numerical simulations as 
shown in Fig. [^ The DOS is symmetric around E = 
and has no singularities, except at the edges E = aE2f3. 




The EDR and DOS of the fundamental vibration v = 1 
of a linear lattice can be obtained analytically by purely 
algebraic methods by making use of certain properties of 
matrices. These properties were first noticed by Nieren- 
berg [31] for lines and further exploited and generalized 
for lines and rings in Refs. [29l|32]. The basic property 
is that the Major ana matrix on a line Ml with matrix 
elements 


= C. (25) 

on the diagonal ”ibs off the main diagonal” has eigenval¬ 
ues given by 


m = Cs2 cos 


skir 


/c = I, 2,..., n. 


(26) 


For nearest-neighbor interactions 5 = 1, the solution can 
be written as 


Ek = a-2/3 cos Ok, Ok 


kiT 


k = ..., n. 


(27) 


This solution applies to both harmonic and anharmonic 
vibrations. The coefficients a and /3 are given by 


FIG. 4. EDR and DOS of the v = 1 vibration of a linear chain 
with i; = I, A = 0.5, A = —I, n = 200, N = 20000. Full red 
lines: Upper panel: Eq. (27), lower panel: Eq. (29). 


2. The overtone vibration v=2 


The spectrum and density of states of the overtone 
vibration can be calculated analytically only in the har¬ 
monic limit. The energies of the v = 2 vibrations are 


Ek^k' =— 2/3 {cosOkcosOk'), k < k\ /c,/c' = I,..., n; 
0- = <“> 

For n large, one can introduce two continuous variables 
^ y = rewrite Eq. (30) as 


E{x, y) = 2a — 2/d(cosx -h cos^) (31) 

with y > X. The DOS can be calculated from Eq. © 
using the method introduced by Bowers and Rosenstock 
[33] and discussed in detail in Sec. IIV A li yielding 


a = —4A 



9{E) 


4K(fcf) 


kl = 4(1 - E)E, 


p = \. 


(28) 


(32) 
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where K is a. complete elliptic integral of the first kind 
and E = {E -2a ^ 4^)/(8/3), 0 < E < 1 . This DOS 
has a logarithmic singularity at E = 1/2. This result is 
confirmed by numerical simulations as shown in Fig. 
where the energy eigenvalues are plotted as a function 




FIG. 5. EDR and DOS of the v = 2 vibration of a linear chain 
in the harmonic limit with A = 0.5, A = —1, n = 100, N = 
20000. Full red lines: upper panel: Eq. (30), lower panel: 
Eq. (I32|. 


of the single index k which labels the eigenvalues in suc¬ 
cessive order, instead of the double indices k, k' used in 
Eq. (|30). 


For anharmonic vibrations, analytic solutions are only 
available for weak ^ 1 and strong ^ ^ 1 couplings 

[29] . Numerical simulations are shown in Fig. The 
EDR exhibits a gap between E^ = 6.4 and Ek = 6 . 6 , 
where the lower sequence contains n eigenvalues and the 
upper one n{n — l)/ 2 , see Eq. (19). The DOS has two 
edge singularities below the gap and one logarithmic sin¬ 
gularity above it. 


IV. SQUARE LATTICE 

For 2 D lattices, the DOS depends on the symmetry 
of the lattice. We consider first the EDR and DOS of a 
square lattice. 

A real-space cell of this lattice is shown in Fig. [^left. 
In this unit cell there are two and only two types of in¬ 
teraction, nearest neighbor and next-to-nearest neighbor. 
For identical atoms, the symmetry of the unit cell, 
imposes the conditions 


41, 

(33) 




FIG. 6. EDR and DOS of the v — 2 vibration of a linear 
anharmonic chain with A = 0.15, A = — 1, n = 100, iV = 10. 
The EDR exhibits a gap between Ek = 6.4 and Ek = 6.6 and, 
accordingly, the DOS vanishes there. 



EIG. 7. The unit cell (left) and supercell (right) of a square 
lattice. 


A real-space supercell of this lattice is shown in Fig.[^ 
right. It consists of four unit cells. The central atom i 
can interact either with its four nearest-neighbors 1, 3, 
5, 7, or with its four next-to-nearest neighbors 2, 4, 6 , 8 . 
The coordinates of the four nearest neighbors are given 

by 


^ = (i,0)> = (0,-i)> = (-1,0), ^ = {0,1), 

a a a a 

(34) 

while those of the next-to-nearest neighbors are 


^ = ( 1 , 1 ), ^ = ( 1 , - 1 ), ^ = (- 1 , - 1 ), ^ = (- 1 , 1 ), 

a a a a 

(35) 

where a is the lattice constant. The symmetry of the 


A^^^ — Ai 2 — A 23 — A 34 — A. 
= Ai3 = A24- 
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lattice imposes the conditions 

= Xii = Xi3 = Xi3 = A^7, 

A^ ^ — Xi 2 — A24 — A26 — A^g* 

The Majorana interaction can then be written as 


M = ^ M'i 

i'f'J) ((hj)) 


ill) 
ij ' 


In the harmonic limit 


M = A(^) ^ (blbj + h.c) + A(") {blbj + h.c.) . 

(38) 

The full Hamiltonian for n atoms in a square lattice is 


H = Eo + AYCi + A'Y C'j + ^ Mlf 


Kj 




A('') ^ mIY- 


(39) 


The operators 


Y Y ( 40 ) 

(hi) ((hi)) 

are also called symmetry adapter operators. 


A. Solutions 

1. The fundamental vibration v=l 

The EDR and the DOS of the fundamental vibration 

= 1 of a square lattice with nearest-neighbor and next- 
to-nearest neighbor interactions and Hamiltonian (38) 
can be obtained analytically by purely algebraic meth¬ 
ods by making use of certain properties of matrices as 
discussed in Sec. IIII AH Eqs. (25) and (26). 

The Hamiltonian matrix is obtained by taking 
matrix elements of H, Eq. (39), in the basis 


|1,0,0,...),|0,1,0,...),-,10,0,...,0,17= 117 ofEq. (jW). 

The operators Q, Ch are diagonal in this basis with 


values given by —4(1 — and —4(^), Eq. (11). The 
interesting operators are and The matrices 

representative of these operators are of the type discussed 
previously and thus can be diagonalized analytically. The 
solution for the fundamental vibration of a square lattice 
was first given by Bowers and Rosenstock [33] in terms of 
functions over an area and subsequently by Nierenberg 
[3T] in terms of functions over a line. In Ref. [33| the 
eigenvalues of and are labelled by two integers 
p, g = 1, 2,..., n and given by 

(n^ \ ^ f PTT \ 

qp, g) = 2 cos-h cos- 

\ n + 1 n + ly 

g) = 4 Tcos cos ^ . (41) 

^ ^ V ^ + 1 n + iy ^ ^ 


Erom these, one can construct the eigenvalues of the Ma¬ 
jorana operator M = X ^^^-h as 


(36) 

m(p, g) = A^^^2 ^ 

(37) 

+ A(^^)4 ( 


- cos - 


piT 

n + 1 

PTT 


n ■ 


1 


cos - 


qiT 

n + 1 
qn 

n + 1 


(42) 


where the minus sign has been introduced to conform 
with the definition in Eq. (13). These solutions can be 


easily verified for the unit cell, n = 4, and for the super¬ 
cell, n = 9. Eor the unit cell, n = 4, the eigenvalues are 
given in the terms of cos(^) and are 

(/) : ± 2 , 0,0 

{II): ±1,±1 (43) 

Eor the supercell, n = 9, they are 

(J) ; ±2\/2, ±V2, ±V2, 0,0,0 
(//) ; ±2, ±2,0,0,0,0,0 (44) 

Eor the unit cell the eigenstates are representations of 
D4/J,. The eigenstates of (I) are the singly degenerate 
representations A, B and the doubly degenerate repre¬ 
sentation E [34] . 

When n is very large, one may replace and 
by continuous variables x and p, respectively, which run 
from 0 to TT, 

m(x, y) = A^^^2 (— cosx — cosp) 

±A^^^)4( — cosxcosp), 0 <x,p< 7r. (45) 

The DOS can be obtained from Eq. ( [4^ . Eor purposes of 
display, it is convenient to consider the shifted function 

m 

f{x, y) = A^'^^2 (2 — cosx — cosy) 

+ A(^^)4(1 — COS X cos p) , 0 < X, p < TT. (46) 

The DOS g{E) is obtained from Eq. (E^ by 


g{E) = T f f 5{E - f{x,y))dxdy. (47) 
Jo Jo 

The DOS is particularly simple when < ^. In 

this case, introducing E^ax = 8A^^^ and E = E/Emax, 
it can be written in terms of an elliptic integral. 


d{E) = 

kl = 


dK{k\) 


(48) 


TT^y(1 + 2\W/\{i)f - 8 (A(^^)/A(^)) E 
4.{l-E)E 

_V_ J_ _ 0 <C -£/ <c 1 

(1 + - 8 E 


Eor nearest-neighbor interactions, A^^^^ = 0, the DOS 
takes the form 


9{E) = 




, ki=A{l-E)E, (49) 
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FIG. 8. Black: DOS of the fundamental vibration u = 1 of 
a square lattice with n = 2500, A = —1,N = 20000. Full 
lines: = 0. Dashed lines: A^^^^/A^^^ = 0.25. Red: 

Eq. ( [4^ . 


FIG. 9. DOS of the fundamental vibration u = 1 of a square 
lattice including the unit cell and the supercell with the num¬ 
ber of sites n as indicated in the insets. Red dashed line: 
Eq. (|49). 


which is symmetric around E = 1/2. The DOS g{E) 
is shown in Fig. for purely nearest-neighbor interac¬ 
tions, A^^^^ = 0, and next-to-nearest neighbor interac¬ 
tions, with A^^^^ = . 

From this ’’reduced” DOS, one can reconstruct the 
’’true” DOS. The energy of the v = 1 states including 
the diagonal terms is given by 


E{p,q) = Eo-4A{l - - 4A'{^) + m{p, q). (50) 


The ’’true” DOS is obtained from Eq. (46) by changing 


^ to and shifting it by the amount Eq — 

4:A (l — — 4^4' (^) — 4A^^^ — 4A^^^\ It is interesting 

to note that the structure of the DOS is a consequence of 
the symmetry D 4 ^h of the unit cell, and is already encoded 
into the DOS of the unit cell and the supercell, as shown 
in Fig. 

This result has been verified by numerical simulations 
as shown in Fig.[^ The DOS of the v = 1 vibration of a 
square lattice has a logarithmic singularity at ^ = 1/2. 

It should be noted that the DOS of the fundamental 
vibration, u = 1, of the 2D lattice is identical to that 
of the overtone vibration u = 2 of the ID lattice. This 
result follows from the mapping of the function f{x^y) 
defined over an area onto a function defined over a line, 
as discussed by Nierenberg m- 


2. The overtone v=2 

While for the fundamental vibration v = 1 the EDR 
and the DOS do not depend on whether or not the crys¬ 
tal is harmonic, for the overtone v = 2 they do. The 
spectrum of states can be calculated by diagonalizing 
the Hamiltonian H in the basis |2,0,0,...), |1,1,0,...), 




EIG. 10. EDR and DOS of the fundamental vibration v — 1 
of a harmonic square lattice with n — 250 0. A — — 1, iV = 
20000, A^^) = 0.5, A^^^^ = 0. Red lines: Eq. (1^). 


|1,0,1,...), ... of Eq. (19), with in total two quanta, in 
terms of the four parameters, A, H', X^^\ X^^^^ and the 
anharmonicity parameter N. Analytic solutions can be 
obtained only in the weak coupling (X/A 1) and strong 
coupling limit (X/A > 1) [2^[32] . 


In order to illustrate the properties of the solutions, 
we consider the case of a highly anharmonic vibration 
with = 10, —4A = 4, —4A' = 0, and only nearest- 
neighbor interactions A^^^ = 0.03, A^^^^ = 0. The DOS 
of u = 2 states has two parts. Both parts diverge as 
n ^ oo, the first part as n and the second as . This is 




















FIG. 11. EDR and DOS of the v = 2 overtone of an anhar- 
monic square lattice with n — 169, A = —1,N = 10, A = 
0.03. 



FIG. 13. The unit cell of the honeycomb lattice, (left), and 
its supercell, (right). 


types of interaction, (I) nearest neighbor, (II) next-to- 
nearest neighbor, and (III) third neighbor. For identical 
atoms the symmetry of the unit cell, imposes the 

condition ([23], page 139) for the interactions, 

= Ai2 = A23 = A34 = A45 = A56 = Ai6, 

A^^^^ = Ai 3 = A24 = A35 = A46 = Ai5 = A26, 

= Ai4 = A 25 = A 36 . (51) 


seen in the numerical simulation in Fig. mi The splitting 
into two pieces containing n and n(n — l)/2 eigenvalues, 
respectively, can already be seen in the unit cell, n = 4, 
and in the supercell, n = 9, as shown in Fig. [T^ 



Ly V _I_^_I__i_ 

-0.9 -0.6 -0.3 0 0.3 

E 


FIG. 12. DOS of the v — 2 overtone of an anharmonic square 
lattice including the unit cell and the supercell with the num¬ 
ber of sites indicated in the insets and A = —1, N = 10, A = 
0.03. 


A real-space supercell of this lattice is shown in Fig.[l^ 
right. It consists of three unit hexagonal cells. The 
central atom i can interact either with its three nearest 
neighbors 1,5,9, or with its six next-to-nearest neighbors 
2,4,6,8,10,12, or with its three third neighbors 3,7,10. 
The coordinates of the nearest neighbors are given by 

^1 = 1 (1, Vs) , ^2 = I (1, - Vs) , 83 = -a ( 1 , 0 ), 

(52) 

while those of the next-to-nearest neighbor, S\ and third 
neighbor, , are given in terms of the lattice vectors 

ai = I (s, Vs) , 02 = I (s, -Vs) , (53) 

as 

= icXi, ^2 ~ ^3 ~ ^ (^2 ~ Otl) (54) 

= <3-1 + ^ 2 ? ^ 2 ” = T ^ 2 ? ^ 3 ” = ~Ci2 + ^ 1 - 

The symmetry of the unit cell of the honeycomb lattice, 
imposes the conditions 

A^^^ = Xii = Xi5 = XiQ, 

X^^^^ = Xi2 = Xi4 = A^6 = A^s = A^io = A^12, 

= Xi3 = Xi7 = XiiQ. (55) 

The Major ana interaction for the honeycomb lattice can 
be written as 


V. HONEYCOMB LATTICE 


A real-space cell of this lattice is shown in the left part 


of Fig. 13 In this unit cell there are three and only three 


{ij) {{ij)) 


{{(ij)}) 


( 56 ) 
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In the harmonic limit, M becomes the Hamiltonian of 
the tight-binding model [3] 

M = ^ (blbj + h.c)j + 

(«d) ((«d)) 

+ ^ {b\bi + h.c). (57) 

((dd))> 

The honeycomb lattice can also be viewed as two inter¬ 
penetrating triangular lattices A and B [21]. In the unit 
cell shown in the left part of Fig.[^ the three atoms 1,5, 
9 belong to A and the three atoms 2, 4, 6 belong to B. Ne¬ 
glecting third neighbor interactions, denoting by a\ and 
bl the boson creation operators at site i on sublattices A 
and B, Eq. ( [4^ can be rewritten as 


M = -t 


E + ^-c-) E (“hi + blbj + h.c?j , 

■ ' {{id)} 


ihj) 


(58) 

where t = and t' = 

The Majorana operator in Eq. (57) is written in terms 
of boson operators. Its structure, however, depends 
only on the symmetry of the lattice. The symme¬ 
try adaptation in Eq. (57) can therefore also be used 
for fermions. Introducing fermion creation operators 
a\ the Hamiltonian has the form 


M = -t + h.c.^ 

— t h.C.^ , (59) 

((hj),cr> 


describing electrons that can hop from one site to the 
other in the interpenetrating triangular lattices A and 
B [3]. The problem of transverse vibrations of a honey¬ 
comb lattice is thus formally identical to that of the band 
structure of the lattice, except for the replacement of bo¬ 
son operators by fermion operators. The only difference 
is that while for bosons one can put any number at each 
site, for fermions one can put only one at each site (or 
two if one takes into accout the spin). 

The full algebraic Hamiltonian for vibrations of the 
honeycomb lattice is 

H = Eo + AY,Ci + A'Y^ Clj + A(^) ^ 

* {ij) 

+ A(") mIY\ (60) 

{{id)) {{{id))) 


There are therefore for this problem three symmetry 
adapter operators 

= ymIP, = Y 

{id) {{id)) 

{{{id))) 


A. Solutions 

1. The fundamental vibration v=l 

The eigenvalues of the symmetry adapter operators 
and for the unit cell are given by 

( J ); ± 2 ,± 1,±1 

{II): ±2,±1,±1 

{III): ±1,±1,±1, (61) 

while for the supercell they are given by 

{!) : ±2.3810, ±1.7320, ±1.5735, ±1,±1, ±0.9246,0 
{II) : 3.6457,3.0861,1,1,0.618,0.428, 

-1, -1, -1, -1.514, -1.618, -1.645, -2 
{III) : ±V2, ±V2, ±1, ±1,0,0.0,0.0. (62) 

The eigenstates are representations of D^h [M]' When 



FIG. 14. The DOS of the v = 1 vibration of the 
lattice, Eq. (64). 


honeycomb 


the number of cells is very large, an analytic solution was 
given by Wallace [21] who solved the equivalent problem 
of fermions on the honeycomb lattice (see also [22]), 

E {kx, ky) = (/c^, ky) + (/c^, ky) (63) 

{kx, ky) = ±xj3Eu{kx,ky) 

^ (^fc: ^y) ^(^fC5 ^1/) 

u{kx, ky) = 2cos27rkya + 4cos7r/cyacos7r/Ca,\/3a, 

where irkxa = x and irkya = y are continuous variables. 
The Hamiltonian used by Wallace differs from M by a 
constant (a shift in E). Hobson and Nierenberg [35] us¬ 
ing a method similar to that described in Sec. IIV A II 
for the square lattice provided an analytic expression 
for the DOS for nearest-neighbor interactions. Intro- 
ducing^ E’max = 6A(^) and E = EjE^^, -1 < .E < 1, 
u = (E + l)/2, the DOS is given in terms of an elliptic 
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integral 


0<cj< 


!<M= (VT^) • 

(v'^) ' 5 ^ 

and the function is symmetric around u; = |, i.e., around 
E = 0, see Fig. M 

The Hamiltonian used by Hobson and Nierenberg dif¬ 
fers from M, and from the Hamiltonian of Wallace by a 
constant. The DOS Eq. (64) should_^ compared with 

The difference 


that of a square lattice, given by Eq. ( |49| ) 
between the two is a consequence of the symmetry of the 
lattice. The major difference between the two DOSs is 
the occurrence of a zero (Dirac zero) at ^ = 0 in the 
honeycomb lattice. The properties of the DOS of the 
honeycomb lattice are due to the symmetry of the 
unit cell and are already encoded into the DOS of the 



-1.2 -0.9 -0.6 -0.3 


0 

E 


0.3 0.6 0.9 1.2 


FIG. 15. DOS of the fundamental vibration u = 1 of the 
honeycomb lattice including the unit cell and the supercell 
with the number of sites n as indicated in the insets. Red 


dashed line: Eq. (64). 


These results have been verified by numerical simula¬ 
tions as shown in Eig. Eor finite n a peak occurs at 
E = 0, which is due to the presence of zigzag edges in 
the lattice structure [36]. 


2. The overtone v=2 


The spectrum (EDR) and density of states (DOS) can 
also be calculated by diagonalizing the Hamiltonian H 
of Eq. (57) in the space v = 2, |2,0,0,...), |1,1,0,...), 
|1, 0,1, 0,...) , ...., in terms of the five parameters H', 
and the anharmonicity parameter, N. 
In order to illustrate the properties of the solutions we 
consider the case of a highly anharmonic vibration with 
N = 10, —4H = 4, —4H' = 0, and only nearest-neighbor 



FIG. 16. EDR and DOS of the fundamental vibration u = 1 
of a honeycomb lattice with n = 4272 sites. 


interactions, = 0.03, = 0. Also here 

the DOS of'T’ = 2 has two parts, as shown in Eig.[^ The 
first part diverges as n and has only edge singularities as 
illustrated in Eig. I^, the second part diverges as n^, and 
has a logarithmic singularity and two shoulders. The 
splitting into these two parts can also be seen in the unit 
cell and in the super cell as shown in Eig. [T^ 




FIG. 17. EDR and DOS of the overtone vibration u = 2 of a 
honeycomb lattice with n = 169, A = — 1, iV = 10, A = 0.03. 
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FIG. 18. Zoom into the small peak in the DOS of the overtone 
vibration i’ = 2 of a honeycomb lattice with n = 169, A = 
— 1, = 10, A = 0.03 shown in Fig. p?| 



V _I_._I_1 "V. _l_ 

-0.9 -0.6 -0.3 0 0.3 

E 


FIG. 19. DOS of the v — 2 overtone of a honeycomb lattice 
including the cell and supercell with the number of sites n 
indicated in the insets and ^ = —1,A^=10, A = 0.03. 


VI. SUMMARY OF DOS AND EDR 


The results of sections m\ |IV| and |V| are summarizec 
in Figs. 20 and 21 For the fundamental vibration u = 1, 
they confirm the conjecture of van Hove m that singu¬ 
larities in the level density occur only in two dimensions. 


VII. ENERGY SURFACES AND PHONON 
DISPERSION RELATIONS 


It is of interest to display explicitly the energy surfaces 
and phonon dispersion relations of infinite-size lattices. 
For the square lattice, the energy surface for nearest- 
neighbor interactions is given by Eq. (46), rewritten, for 


FIG. 20. Summary of the EDRs of transverse crystal vibra¬ 
tions. Left: i’ = 1, right: v — 2. 



FIG. 21. Same as Fig.j^for the DOEs. 


2A^^^ = 1, as 

E{kx, ky) = (2 — cosirkxa — cosnkya). (65) 


This energy surface is shown in Eig.[^ The correspond¬ 
ing phonon dispersion relation along the boundary of the 
first irreducible Brillouin zone is shown in Eig. The 
special points F(0,0), M (^,0), X (^, ^) are indicated 
in the density plot, shown in the right part of Eig. 

Eor the hexagonal lattice, the energy surface for 
nearest-neighbor interaction is given by Eq.(62) rewritten 
as [3] 


E{kx, ky) = ±^3 + f{kx,ky) (66) 

f{kx,ky) = 2cos (^VSkyO^ +4cos cos ■ 

This energy surface is shown in Eig. Because of the 
± sign in Eq.(66), it consists of two sheets which touch 
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FIG. 22. The left panel shows the energy surface of the v — 1 
vibrations, Eq. (65) of an infinite-size square lattice, where 
the lattice constant a was set to unity. The right panel shows 
the corresponding density plot in the quasimomentum plane 
(/Ccc, ky) with the isofrequency lines shown as dark lines. 



FIG. 23. Gomputed phonon dispersion relation of an infinite- 
size square lattice. 




FIG. 24. The left panel shows the two sheets of the energy 
surface of the v — 1 vibrations, Eq. (66) of an infinite-size 
hexagonal lattice, where a was set to unity. They touch each 
other conically at the corners of the first Brillouin zone. The 
right panel shows the corresponding density plot in the quasi- 
momentum plane {kx,ky) with the isofrequency lines shown 
as dark lines. 


conically at the corners of the first Brillouin zone. This 


± sign is a consequence of the fact that the hexagonal 
lattice can be viewed as two interpenetrating triangu¬ 
lar lattices as discussed in Sec. m This is the crucial 
property that makes the hexagonal lattice so different 
from other lattices. The phonon dispersion relation is 
shown in Fig. 25 The special points F(0,0), M 



FIG. 25. Gomputed phonon dispersion relation of an infinite- 
size hexagonal lattice. 


K (l^, are indicated in the density plot, shown in 

the right part of Fig. In Figs. and one can see 
clearly the symmetry of the unit cells (square) and 
D^h (hexagonal). 

As mentioned in Sec. m the structure of the algebraic 
Hamiltonian applies to both bosons and fermions, and 
it depends only on the symmetry of the lattice. Thus, 
the energy surfaces and dispersion relations, Eqs. (65) 
and (66), apply also to electrons in a square and hexag¬ 
onal lattice. In this case, k = {kx, ky) represents 
the quasimomentum of electrons. This is what makes 
graphene so special. 


VIII. COMPARISON WITH DATA IN 
MICROWAVE PHOTONIC CRYSTALS 

Superconducting microwave resonators have been used 
for two decades as analog systems for the study of 
quantum phenomena in high resolution measurements 
El 1 0 |38]. Photonic crystals are the optical analog of 
a solid. Both concepts can be combined into ’’microwave 
photonic crystals” which offer the opportunity to perform 
high precision measurements of the excitation spectrum, 
and thus to study the FDR and DOS of solids. In partic¬ 
ular, recently two-dimensional hexagonal structures were 
realized by squeezing a photonic crystal, which was com¬ 
posed of several hundreds of metallic cylinders forming a 
triangular lattice, between two metal plates [201 EH So]. 

In order to investigate the FDR and the DOS of hexag¬ 
onal lattices of varying shapes high-precision experiments 
were performed with flat, superconducting microwave 
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resonators with the forms of a rectangle and the African 
continent [41], respectively. Photographs are shown in 
the left parts of Figs. [^ and ^7\ Both resonators con- 



FIG. 26. (Color online) Left part: Photograph of the rectan¬ 
gular microwave photonic crystal, which contains 888 metal 
cylinders arranged on a triangular grid. The top plate was 
shifted with respect to the bottom one for presentational rea¬ 
sons. Right part: Schematic view of the triangular lattice. 
The red (gray) and blue (dark gray) dots mark the voids be¬ 
tween the cylinders that form the hexagonal lattice. Adopted 
from [38] 



FIG. 27. (Color online) Same as in Fig. 26 but for an Africa 
billiard. The top plate was removed. Adopted from [38] 


sisted of a basin and a lid, that were made from brass 
plates and then lead plated to attain superconductivity 
at the liquid helium temperature of 4.2 K. Note that 
the critical temperature of lead is Tc = 7.2 K. For 
the construction of the photonic crystals located inside 
the basins ~ 900 cylinders were milled out of the bot¬ 
tom plate. The cylinders were arranged on a triangu¬ 
lar grid, as indicated in the schematic views in Figs. ^6\ 
and[27l Then the voids at the centers of the cells formed 
by, respectively, three of them yield a hexagonal config¬ 
uration. The height of the resonators was d = 3 mm 
and the range of excitation frequences / of the mi¬ 
crowaves that were coupled into the resonator was chosen 
as 0 < / < fmax = c/2d so the electric field vector was 
perpendicular to the top and bottom plates of the res¬ 
onators. Then the microwaves inside the resonators are 
governed by the scalar Helmholtz equation for the electric 
field strength with Dirichlet boundary conditions at the 
walls of the cylinders and the basin and the microwave 
photonic crystals correspond to experimental realizations 
of hexagonal lattices. Consequently, the model Eq. (59) 
should be applicable. 

We determined altogether 1651 and 1823 resonance fre¬ 
quencies in the first two bands of the rectangular and 
the Africa-shaped microwave photonic crystal, respec¬ 


tively. The EDRs, fk vs. k, and the DOSs p (/) are 

shown in Eigs. [^ and [^ A fit to the DOS with the 

algebraic Hamiltonian Eq. (60) with nearest-neighbor 
coupling, and second- and third-neighbor cou¬ 
plings, and yielded the red dashed curves 

in Eig. [^ The best fit values are for the rectangular 
crystal = 4.573 GHz, = -0.284 GHz, = 

0.104 GHz, and for the crystal with the shape of Africa 

A^^) = 2.198 GHz, A^^^) = 0.107 GHz, = -0.013. 

The data exhibit both the van Hove singularities and 
a vanishing DOS at the Dirac frequency [9]. 



FIG. 28. Experimental EDR for the rectangular photonic 
crystal with 1656 sites (upper panel) and the Africa-shaped 
one with 1823 sites (lower panel). 


The DOS diverges logarithmically at the van Hove sin¬ 
gularities only for 2-dimensional, periodic structures of 
infinite extent. In the crystals used in the experiments, 
however, the sharp peaks in the DOS have a finite height 
^max determined for the experimental DOS of 
the two microwave photonic crystals and a third, smaller 
one which contained only 267 cylinders, and also per¬ 
formed numerical studies for photonic crystals of various 
sizes with the shapes of a rectangle and of Africa. Eor 
a comparison of these results we rescaled the frequencies 
such that the distance between the van Hove singularities 
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T-1-1-1-1-1-1-1-1-1-1-1-1-1-r 



/(GHz) 



/(GHz) 


FIG. 29. Same as Fig. for the experimental DOS. The 
red dashed lines show the hts of the DOS deduced from the 
Hamiltonian in Eq. (60) to the experimental ones. 


was 2 for all systems. The experimental and numerical 
studies revealed that the maxima of the DOS, be¬ 

have like 


^max ^ [\og{Nc) + b] (67) 

with Nc denoting the number of unit cells, i.e., of 
hexagons formed by the voids in the photonic crystal. 
Here, a and b are fit parameters, where the former takes 
a similar value a ^ 0.141 — 0.155 ~ for all cases, i.e., 
it seems to be universal. 


IX. CONCLUSIONS 

In this article. Id (transverse out-of-plane) z-vibrations 
have been studied. For the analysis of the experimen¬ 
tal data, one needs to generalize the method to 2 d (lon¬ 
gitudinal in-plane) xp-vibrations. In the algebraic ap¬ 
proach, this is done by introducing the Lie algebra g = 
u{3) [42l[43]. The combined algebra at each site is then 
1^(2) 0 u{3). The application of the algebraic method to 


0.5 


0.4 


CL 

0.2 


2.5 3 3.5 4 

log(A^,) 

FIG. 30. The height of the van Hove peaks in the DOS 
of photonic crystals of varying size versus the logarithm of the 
number of hexagonal cells Nc. Black dots: Numerical results 
for photonic crystals with the shape of Africa. The dashed 
line shows a straight-line ht a log (Nc) + b through these data, 
where the resulting slope equaled a = 0.141. Red squares: 
Numerical results for photonic crystals with the shape of a 
rectangle. Turquoise triangles: Experimental results. 



longitudinal vibrations and to combinations of transverse 
and longitudinal vibrations will be reported in a subse¬ 
quent publication. Also, for those situations in which the 
lattice is not composed of identical units X — X — X —..., 
but it has alternating units X — Y — X — Y — ... ^ one 
needs to compute the phonon dispersion relation both 
for the optical and acoustic branches. This is easily done 
within the framework of the algebraic theory discussed 
here. 


Finally, we have presented a new method for calcu¬ 
lating the energy spectrum (EDR) and density of states 
(DOS) of vibrations of solids, both harmonic and anhar- 
monic, and applied it to the study of the EDRs and the 
DOSs of ID linear chains, and 2 D square and hexago¬ 
nal lattices. Our results have been compared with data 
obtained in microwave photonic crystals for 2D hexago¬ 
nal lattices. These data show the expected occurrence of 
both van Hove singularities and Dirac zeros. The method 
can be easily extended to other 2D-lattices, with symme¬ 
try of the unit cell other than D 4 ^h and and in fact 
to any two-dimensional structure planar and non-planar, 
as for example fullerene, Ceo, with icosahedral symmetry 
J/i. The algebraic method can also be used to calculate 
the response of a solid to infrared (IR) and Raman (R) 
radiation. For ID systems, this response was studied in 
[29] . For 2D systems it remains to be done. 


The method is also well suited to study the general 
Hubbard boson Hamiltonian of Eq. (24). Because of 
the formal equivalence between the boson Hamiltonian, 
Eq. (58), and the fermion Hamiltonian, Eq. (59), it can 
also be used to study the band structure in ID, 2D square 
and hexagonal lattices in the tight-binding model. The 











DOS of the fundamental vibration, u = 1, of an hexag¬ 
onal lattice is identical to the DOS of electrons in the 
same lattice. 
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